!!****m* ABINIT/m_paw_slater
!! NAME
!!  m_paw_slater
!!
!! FUNCTION
!!  This module defines objects and procedures to evaluate Slater-like integrals
!!  using real spherical Harmonics.
!!
!! COPYRIGHT
!! Copyright (C) 2008-2020 ABINIT group (MG)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!!
!! NOTES
!!  * Routines tagged with "@type_name" are tightly connected to the definition of the data type.
!!    Tightly connected means that the proper functioning of the implementation relies on the
!!    assumption that the tagged procedure is consistent with the type declaration.
!!    Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
!!    that all the tightly connected routines are changed accordingly to accommodate the modification of the data type.
!!    Typical examples of tightly connected routines are creation, destruction or reset methods.
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

MODULE m_paw_slater

 use defs_basis
 use m_abicore
 use m_errors
 use m_splines

 use m_fstrings,     only : basename
 use m_paw_atomorb,  only : atomorb_type, init_atomorb, print_atomorb, destroy_atomorb, get_overlap
 use m_crystal,      only : crystal_t
 use m_paw_io,       only : pawio_print_ij
 use m_pawang,       only : pawang_type
 use m_paw_sphharm,   only : realgaunt
 use m_pawrad,       only : pawrad_type, pawrad_free, pawrad_isame, &
&                           pawrad_deducer0, simp_gen, calc_slatradl
 use m_pawtab,       only : pawtab_type
 use m_pawrhoij,     only : pawrhoij_type
 use m_paw_lmn,      only : make_kln2ln, make_klm2lm, make_indln, klmn2ijlmn

 implicit none

 private

 public :: paw_sigxcore            ! The onsite matrix elements of the Fock operator generated by (closed) core shells.
 public :: paw_mkdijexc_core       ! Calculate the onsite matrix element of the Fock operator generated by the core.
 public :: paw_dijhf               ! Compute the onsite D_{ij} strengths of the exchange parth of the self energy.
!!***

!!****t* m_paw_slater/slatrad_t
!! NAME
!!  slatrad_t
!!
!! FUNCTION
!!  Object used to store radial integrals of the form.
!!
!!  $ F_{ijkl}^L = \dfrac{4\pi}{2L+1} \int u_i(1) u_j(1) \dfrac{r_<^L}{r_>^{L+1}} u_k(2) u_l(2) d1d2 $
!!
!!  for a given quadruple (i,j,k,l) as a function L \in [L_min, L_max].
!!  i,j,k,l are shorthand indeces for (nn,ll) quantum numbers.
!!
!! NOTES
!!   Basic symmetry properties:
!!   1) invariant under the exchange (i<-->j) and (k<-->l).
!!   2) invariant under the exchange (i,j) <--> (k,l).
!!
!!  Memory saving is achieved by storing the upper triangle of the (ij) (kl) matrix
!!  and, for each dimension, only the upper triangle of the two matrices (iln,jln) (kln,lln).
!!
!!  Some matrix elements will never contribute to <ij|1/|1-2||kl> due to selection rules
!!  introduced by the integration of the angular part.
!!
!! SOURCE

 type,public :: slatrad_t

  integer :: iln,jln,kln,lln
  ! The (l,n) indeces associated to the partial waves.

  integer :: lslat_min
  ! Min l+1 in the expansion of the Coulomb potential.

  integer :: lslat_max
  ! Max l+1 in the expansion of the Coulomb potential.

  integer :: nintgl
  ! The number of non-zero integrals stored in intgl.

  integer,allocatable :: intgl_select(:)
  ! intgl_select(lslat_min:lslat_max)
  ! Index of the non-zero integrals in intgl, 0 if intgl has not been
  ! calculated thanks to selection rules coming from the angular integration.

  real(dp),allocatable :: intgl(:)
  ! intgl(1:nintgl)
  ! The integrals:
  ! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) \phi_{\nj\lj}(1) \dfrac{r_<^L}{r_>^{L+1}} \phi_{\nk\lk}(2) \phi_{\nl\ll} d1d2
  ! for given (i,j,k,l) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).

 end type slatrad_t

 public :: slatrad_init      ! Creation method
 public :: slatrad_free      ! Free memory
!!***

 interface slatrad_free
   module procedure slatrad_free_0D
   module procedure slatrad_free_1D
 end interface slatrad_free

!----------------------------------------------------------------------

!!****t* m_paw_slater/slatang_cshell_t
!! NAME
!!  slatang_cshell_t
!!
!! FUNCTION
!!  Object used to store:
!!   $ F^{lsl,lc}_{li,lj,mi,mj} = sum_{msl mc} \<li mi|lsl msl;lc mc\> \<lsl msl;lc mc| lj mj\> $
!!  This (less general) type of radial integral is needed to evaluate the Exchange term generated
!!  by a closed-shell atom. In the equation, (lc,mc) are the set of angular quantum number associated
!!  to (closed) core electrons while  (lsl,msl) comes from the expansion of 1/|r1-r2|.
!!  Since the F is invariant under exchange of i and j we use an array of structures indexed
!!  by kln = (iln,jln) in packed form.
!!
!! SOURCE

 type, public :: slatang_cshell_t

  integer :: nsggsel
  ! Number of non null matrix elements

  integer :: lslat_max
  ! Max l+1 in the expansion of the Coulomb potential

  integer :: lslat_min
  ! Min l+1 in the expansion of the Coulomb potential

  integer :: lc_max
  ! Max l+1 for orbitals summed over (usually core orbitals)

  integer,allocatable :: sggselect(:,:)
  ! sggselect(lslat_max,lc_max)
  ! Index of non null sgg, 0 if sgg is zero by symmetry.

  real(dp),allocatable :: sgg(:)
  ! sgg(nsggsel)
  ! Non null matrix elements in packed form. The index is given by sggselect.

 end type slatang_cshell_t

 public :: slatang_cshell_init  ! Creation method for slatang_cshell_t.
 public :: slatang_cshell_free  ! Destruction method for the slatang_cshell_t.
!!***

!----------------------------------------------------------------------

!!****t* m_paw_slater/slatrad_cshell_t
!! NAME
!!  slatrad_cshell_t
!!
!! FUNCTION
!!  Object used to store the set of radial integrals:
!!
!!  $ \dfrac{4\pi}{2L+1} \times
!!    \int \phi_{\ni\li}(1) \phic_{\nc\lc}(1) \dfrac{r_<^L}{r_>^{L+1}} \phic_{\nc\lc}(2) \phi_{\nj\lj} d1d2 $
!!
!!  for given (in,il) and (jn,jl) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
!!  This (less general) type of radial integral is needed to evaluate the Exchange term generated
!!  by a closed-shell atom. In the equation, (\nc,\lc) are the set of angular quantum number associated
!!  to core electrons while  (lsl,msl) comes from the expansion of 1/|r1-r2|.
!!  Since the F is invariant under exchange of i and j we use an array of structures indexed
!!  by kln = (iln,jln) in packed form.
!!
!! SOURCE

 type,public :: slatrad_cshell_t

  integer :: lnc_size
  ! Number of (n,l) channel for core orbitals.

  integer :: lslat_max
  ! Max l+1 in the expansion of the Coulomb potential

  integer :: lslat_min
  ! Min l+1 in the expansion of the Coulomb potential

  integer :: nrlphic_int
  ! The number of non-zero integrals stored in rlphic_int.

  integer,allocatable :: rlphic_select(:,:)
  ! rlphic_select(lslat_max,lnc_size)  TODO should be allocated with lslat_min:lslat_max
  ! Index of the non-zero integrals in rlphic_int, 0 if rlphic_int has not been
  ! calculated thanks to selection rules coming from the angular integration.

  real(dp),allocatable :: rlphic_int(:)
  ! rlphic_int(1:nrlphic_int)
  ! The integrals:
  ! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) \phic_{nc\lc}(1) \dfrac{r_<^L}{r_>^{L+1}} \phic_{nc\lc}(2) \phi_{\nj\lj} d1d2
  ! for given (in,il) and (jn,jl) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).

 end type slatrad_cshell_t


CONTAINS  !========================================================================================
!!***

!!****f* m_paw_slater/slatang_cshell_init
!! NAME
!!  slatang_cshell_init
!!
!! FUNCTION
!!  Initialize the structure slatrad_cshell_t containing radial integrals, see below.
!!
!! INPUTS
!!  l_max= max. value of ang. momentum l+1;
!!   Gaunt coeffs up to [(2*l_max-1,m),(l_max,m),(l_max,m)] are computed
!!  lc_max=Max Lc+1 for core states used to contrusct \Sigma_x^\core.
!!  lm2_size=Number of symmetrix elements in the (l,m) basis set.
!!   NB: lm2_size = (l_max**2)*(l_max**2+1)/2.
!!  klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm and jlm for each klm=(ilm,jlm)
!!  where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection.
!!
!! OUTPUT
!!  Slatang3l(lm2_size) <type(slatang_cshell_t)> = Object storing :
!!
!!   $ F^{ilsl,ilc}_{klm} = sum_{msl,mc}  <li mi|lsl msl;lc mc> <lsl msl;lc mc|lj mj> $
!!
!!  where klm = runs over the upper triangle of the ((il,im),(jl,jm)) matrix.
!!  ilc runs from 1 up to lc_max and |li-lc| <= lsl <= |li+lc|
!!
!! NOTES
!!  Selection rules for F
!!   1) mi = mj
!!   2) In the case of closed shells, one sums for all possible mc"s values from -lc up to +lc.
!!      In this particular case, one can use the symmetry properties of Clebsch-Gordan
!!      coefficients to show that F is non null only if li==lj. In particular, F can be rewritten as:
!!
!!      $ F^{ilsl,ilc}_{klm} = \delta{li,lj}\delta{mi,mj} \times
!!         \Gaunt^{lsl,0}_{lc,0;li,0} \sqrt{ \dfrac{(2*lc+1) (2*lsl+1)}{4\pi*(2*li+1)} } $
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: l_max,lc_max,lm2_size
!arrays
 integer,intent(in) :: klm2lm(6,lm2_size)
 type(slatang_cshell_t),intent(out) :: Slatang3l(lm2_size)

!Local variables-------------------------------
!scalars
 integer :: ilm,ilm0,jlm,lgnt_max,ngnt,ilsl,ilc,lc,ilm0c,lsl
 integer :: klm_ci,ilm0sl,li,il,jl,im,jm,ig000
 integer :: klm,k0lm_i,k0lm_j,k0lm_c,nsggsel,lslat_max,lslat_min
 real(dp) :: dum
!arrays
 integer,allocatable :: gntselect(:,:)
 real(dp),allocatable :: realgnt(:),tmp_sgg(:)
 real(dp),allocatable :: g000(:,:,:)

! *************************************************************************

 !@slatang_cshell_t

 ! * Calculate $\Gaunt^{lsl,msl}_{lc,mc;li,mi}$
 lgnt_max = MAX(l_max,lc_max)
 ABI_MALLOC(  realgnt,((2*lgnt_max-1)**2*(lgnt_max)**4))
 ABI_MALLOC(gntselect,((2*lgnt_max-1)**2, lgnt_max**2*(lgnt_max**2+1)/2))

 call realgaunt(lgnt_max,ngnt,gntselect,realgnt)

 ! Below we need $\Gaunt_{lsl,0}_{lc,0;li,0}$
 ABI_MALLOC(g000,(2*lgnt_max-1,lc_max,l_max))
 g000 = zero

 do ilsl=1,2*lgnt_max-1
   lsl    = ilsl-1
   ilm0sl = 1+lsl**2+lsl
   do il=1,l_max
     li    = il-1
     ilm0  = 1+li**2+li
     k0lm_i = ilm0 *(ilm0-1)/2
     do ilc=1,lc_max
       lc    = ilc-1
       ilm0c = 1+lc**2+lc
       k0lm_c= ilm0c * (ilm0c-1)/2
       if (ilm0c > ilm0) then
         klm_ci = k0lm_c + ilm0
       else
         klm_ci = k0lm_i + ilm0c
       end if
       ig000 = gntselect(ilm0sl,klm_ci) ! Index of $\Gaunt_{lsl,0}_{lc,0;li,0}$
       if (ig000 > 0) g000(ilsl,ilc,il)=realgnt(ig000)
     end do
   end do
 end do

 ABI_FREE(realgnt)
 ABI_FREE(gntselect)

 ! === Loop over klm channels in packed form ===
 do klm=1,lm2_size
   il = klm2lm(1,klm); im = klm2lm(3,klm)
   jl = klm2lm(2,klm); jm = klm2lm(4,klm)

   nsggsel=0
   lslat_min = 1 !FIXME find better way
   lslat_max = il+lc_max-1

   Slatang3l(klm)%lslat_min = lslat_min
   Slatang3l(klm)%lslat_max = lslat_max
   Slatang3l(klm)%lc_max    = lc_max

   ABI_MALLOC(tmp_sgg,(lslat_max*lc_max))
   tmp_sgg = zero

   ! === Calculate F^{lsl,lc}_{li,mi;lj,mj} ===
   ! * Selection rule: mi = mj and li==lj
   if (im == jm .and. il==jl) then
     li  = il-1
     ilm = klm2lm(5,klm); k0lm_i = ilm *(ilm -1)/2
     jlm = klm2lm(6,klm); k0lm_j = jlm *(jlm -1)/2

     ABI_MALLOC(Slatang3l(klm)%sggselect,(lslat_max,lc_max))
     Slatang3l(klm)%sggselect = 0

     do ilsl=lslat_min,lslat_max
     !% do ilsl=lslat_min,lslat_max,2
       lsl = ilsl-1
       do ilc=1,lc_max
         lc = ilc-1
         dum = SQRT( (two*lc+1)*(two*lsl+1) / (four_pi*(two*li+1)) ) * g000(ilsl,ilc,il)
         if (ABS(dum)>=tol12) then ! * Store results and progressive index if non null.
           nsggsel = nsggsel + 1
           tmp_sgg(nsggsel) = dum
           Slatang3l(klm)%sggselect(ilsl,ilc) = nsggsel
         end if
       end do !ilc
     end do !ilsl
   end if ! Selection rule li=lj and mi == mj
   !
   ! * Finalize the object.
   Slatang3l(klm)%nsggsel = nsggsel
   if (nsggsel > 0) then
     ABI_MALLOC(Slatang3l(klm)%sgg,(nsggsel))
     Slatang3l(klm)%sgg = tmp_sgg(1:nsggsel)
   end if
   ABI_FREE(tmp_sgg)
 end do !klm

 ABI_FREE(g000)

#if 0
! Debugging code
 do klm=1,lm2_size
   if (Slatang3l(klm)%nsggsel>0) then
     il  = klm2lm(1,klm)
     jl  = klm2lm(2,klm)
     im  = klm2lm(3,klm)
     jm  = klm2lm(4,klm)
     write(std_out,*)"--for li, mi",il-1,im-il
     lslat_min = Slatang3l(klm)%lslat_min
     lslat_max = Slatang3l(klm)%lslat_max

     do ilc=1,lc_max
       do ilsl=lslat_min,lslat_max
         ii = Slatang3l(klm)%sggselect(ilsl,ilc)
         if (ii>0) write(std_out,*)"   lc, lslat, sgg",ilc-1,ilsl-1,Slatang3l(klm)%sgg(ii)
       end do
     end do
   end if
 end do
#endif

end subroutine slatang_cshell_init
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/slatang_cshell_free
!! NAME
!!  slatang_cshell_free
!!
!! FUNCTION
!!  Free the dynamic memory allocated in a structure of type slatang_cshell_t
!!
!! SIDE EFFECTS
!!  Slatang3l(lm2_size) <type(slatang_cshell_t)> = Object containing radial integrals
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatang_cshell_free(Slatang3l)

 implicit none

!Arguments ------------------------------------
!scalars
 type(slatang_cshell_t),intent(inout) :: Slatang3l(:)

!Local variables-------------------------------
 integer :: ii
! *********************************************************************

 !@slatang_cshell_t
 do ii=1,SIZE(Slatang3l)
   if (allocated(Slatang3l(ii)%sggselect)) then
     ABI_FREE(Slatang3l(ii)%sggselect)
   end if
   if (allocated(Slatang3l(ii)%sgg)) then
     ABI_FREE(Slatang3l(ii)%sgg)
   end if
 end do

end subroutine slatang_cshell_free
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/slatrad_cshell_free
!! NAME
!!  slatrad_cshell_free
!!
!! FUNCTION
!!  Free the dynamic memory allocated in a structure of type slatrad_cshell_t
!!
!! SIDE EFFECTS
!!  Slatrad3l(ln2_size) <type(slarad3l_type)> = Object containing radial integrals
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatrad_cshell_free(Slatrad3l)

 implicit none

!Arguments ------------------------------------
!scalars
 type(slatrad_cshell_t),intent(inout) :: Slatrad3l(:)

!Local variables-------------------------------
 integer :: ii
! *********************************************************************

 !@slatrad_cshell_t
 do ii=1,SIZE(Slatrad3l)
   if (allocated(Slatrad3l(ii)%rlphic_select)) then
     ABI_FREE(Slatrad3l(ii)%rlphic_select)
   end if
   if (allocated(Slatrad3l(ii)%rlphic_int)) then
     ABI_FREE(Slatrad3l(ii)%rlphic_int)
   end if
 end do

end subroutine slatrad_cshell_free
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/slatrad_cshell_init
!! NAME
!!  slatrad_cshell_init
!!
!! FUNCTION
!!  Initialize the structure storing the radial part of Slater"s integrals.
!!
!! INPUTS
!!  ln2_size=Number of symmetrical (l,n) channels
!!  Pawrad<pawrad_type>=paw radial mesh and related data
!!  Pawtab<pawtab_type>=paw tabulated starting data
!!  Atm<atomorb_type>=Object containing core orbitals.
!!  Atmrad<pawrad_type>=paw radial mesh and related data for the atom.
!!  kln_mask
!!
!! OUTPUT
!!  Slatrad3l<slatrad_cshell_t>=The object completely initialized.
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad,kln_mask)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: ln2_size
!arrays
 integer,optional,intent(in) :: kln_mask(ln2_size)
 type(atomorb_type),intent(in) :: Atm
 type(pawrad_type),target,intent(in) :: Atmrad,Pawrad
 type(pawtab_type),target,intent(in) :: Pawtab
 type(slatrad_cshell_t),intent(out) :: Slatrad3l(ln2_size)

!Local variables ---------------------------------------
!scalars
 integer :: cmesh_size,dmesh_size
 integer :: il,iln,ilnc,isl,in,jl,jln,jn,kln,ll,lnc_size
 integer :: lslat_max,lslat_min,lc_max,nintg
 integer :: lmn_size,lmn2_size,do_spline,ln_size,whichdenser,isppol
 real(dp) :: intg,intg1,ybcbeg,ybcend
 logical :: hasameq
!arrays
 integer,allocatable :: kln2ln(:,:)
 integer, pointer :: indklmn(:,:),indlmn(:,:)
 real(dp),allocatable :: ff1(:),ff2(:),tmp_integrals(:)
 real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:)
 real(dp),allocatable,target :: phi_spl(:,:)
 real(dp),allocatable :: der(:),ypp(:)
 real(dp),ABI_CONTIGUOUS pointer :: crad(:),drad(:),phi_in(:)

! *************************************************************************

 ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
 if (PRESENT(kln_mask)) then
   MSG_ERROR("kln_mask is present")
 end if

 !@slatrad_cshell_t
 lmn_size  = Pawtab%lmn_size
 lmn2_size = Pawtab%lmn2_size
 ln_size   = Pawtab%basis_size

 lnc_size = Atm%ln_size
 lc_max   = Atm%l_max

 call pawrad_isame(Atmrad,Pawrad,hasameq,whichdenser)

 do_spline=0
 if (.not.hasameq) then
   do_spline=1
   if (whichdenser/=1) &
&    MSG_COMMENT("Pawrad is denser than Atmrad!")
 else
   ABI_CHECK(whichdenser==1,"Pawrad is denser than Atmrad!")
 end if

 dmesh_size = Atmrad%mesh_size
 cmesh_size = Pawtab%mesh_size

 drad  => Atmrad%rad(1:dmesh_size)
 crad  => Pawrad%rad(1:cmesh_size)

 ! === Spline valence basis set onto core mesh (natural spline) ===
 if (do_spline==1) then
   MSG_COMMENT("Splining in init_slatrad3l")
   ABI_MALLOC(phi_spl,(dmesh_size,ln_size))
   ABI_MALLOC(der,(cmesh_size))
   ABI_MALLOC(ypp,(cmesh_size))

   do iln=1,ln_size
     phi_in => Pawtab%phi(:,iln)
     ypp(:) = zero; ybcbeg = zero; ybcend = zero
     call spline(crad,phi_in,cmesh_size,ybcbeg,ybcend,ypp)
     call splint(cmesh_size,crad,phi_in,ypp,dmesh_size,drad,phi_spl(:,iln))
   end do

   ABI_FREE(der)
   ABI_FREE(ypp)
 end if

 indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
 indklmn => Pawtab%indklmn(1:8,1:lmn2_size)

 ABI_MALLOC(kln2ln,(6,ln2_size))

 call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)

 ABI_MALLOC(ff1,(dmesh_size))
 ABI_MALLOC(ff2,(dmesh_size))

 ! * Loop over the upper triangle of the [(in,il), (jn,il)] matrix.
 ABI_CHECK(Atm%nsppol==1,"nsppol==2 not tested")

 do isppol=1,Atm%nsppol
   do kln=1,ln2_size
     il  = kln2ln(1,kln)
     jl  = kln2ln(2,kln)
     in  = kln2ln(3,kln)
     jn  = kln2ln(4,kln)
     iln = kln2ln(5,kln)
     jln = kln2ln(6,kln)

     lslat_max = MAX((il+lc_max),(jl+lc_max))       - 1   ! These are indeces, not l-values.
     !lslat_min = MIN(ABS(il-lc_max),ABS(jl-lc_max)) + 1
     lslat_min = 1 ! FIXME find better way

     Slatrad3l(kln)%lnc_size    = lnc_size
     Slatrad3l(kln)%lslat_min   = lslat_min
     Slatrad3l(kln)%lslat_max   = lslat_max

     Slatrad3l(kln)%nrlphic_int = 0

     ABI_MALLOC(Slatrad3l(kln)%rlphic_select,(lslat_max,lnc_size))
     Slatrad3l(kln)%rlphic_select(:,:) = 0

     !if (PRESENT(kln_mask)) then  !FIXME THIS IS WRONG, move it below in case
     ! if (kln_mask(kln)==0) CYCLE
     !end if

     if (do_spline==1) then
       MSG_COMMENT("Performing spline of valence phi")
       phi_i => phi_spl(:,iln)
       phi_j => phi_spl(:,jln)
     else
       phi_i => Pawtab%phi(:,iln)
       phi_j => Pawtab%phi(:,jln)
     end if

     ! * Loop over (n,l) channels for Atom orbitals
     ABI_MALLOC(tmp_integrals,(lslat_max*lnc_size))
     tmp_integrals(:) = zero
     nintg=0

     do ilnc=1,lnc_size
       ! phicore => Atm%phi(:,ilnc,isppol)
       ff1 = phi_i * Atm%phi(:,ilnc,isppol)
       ff2 = phi_j * Atm%phi(:,ilnc,isppol)
       do isl=lslat_min,lslat_max ! L coming from Coulomb expansion
         ll = isl-1
         call calc_slatradl(ll,dmesh_size,ff2,ff1,Atmrad,intg1)
         call calc_slatradl(ll,dmesh_size,ff1,ff2,Atmrad,intg)

         !call calc_slatradl(ll,cmesh_size,ff2,ff1,Pawrad,intg1)
         !call calc_slatradl(ll,cmesh_size,ff1,ff2,Pawrad,intg)

         if (ABS(intg1-intg)>tol6) write(std_out,*)"DEBUG ",ll,il,in,jl,jn,intg1,intg

         ! * Store results
         if (ABS(intg)>=tol12) then
           nintg = nintg +1
           Slatrad3l(kln)%rlphic_select(isl,ilnc) = nintg
           tmp_integrals(nintg) = intg
         end if
       end do !ll
     end do ! ilnc

     ! Finalize the object
     Slatrad3l(kln)%nrlphic_int = nintg
     ABI_MALLOC(Slatrad3l(kln)%rlphic_int,(nintg))
     if (nintg>0) Slatrad3l(kln)%rlphic_int(:) = tmp_integrals(1:nintg)

     ABI_FREE(tmp_integrals)
   end do !kln
 end do !isppol

 ABI_FREE(ff1)
 ABI_FREE(ff2)
 ABI_FREE(kln2ln)

 if (do_spline==1)  then
   ABI_FREE(phi_spl)
 end if

end subroutine slatrad_cshell_init
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/paw_sigxcore
!! NAME
!!  paw_sigxcore
!!
!! FUNCTION
!!  Calculate the integrals:
!!  \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) orb_{nl}(1) \dfrac{r_<^L}{r_>^{L+1}} orb_{nl}(2) \phi_{\nj\lj} d1d2!!
!!  for given (in,il) and (jn,jl) as a function of (nc,lc) and L = |il-jl|, |il-jl|+2, ..., |il+il|
!!
!! INPUTS
!!  cplex_dij=1 if dijexc_core is real, 2 if they are complex
!!  lmn2_size=Number of (klmn) channels
!!  ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
!!  Pawtab<pawtab_type>=paw tabulated starting data
!!  Atm<atomorb_type>=Structure containing core orbitals
!!  Atmrad<pawrad_type>=The radial mesh for core orbitals
!!
!! OUTPUT
!!  dijexc_core(cplex_dij*lmn2_size,ndij)
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine paw_sigxcore(cplex_dij,lmn2_size,ndij,Pawrad,Pawtab,Atm,Atmrad,dijexc_core)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: lmn2_size,cplex_dij,ndij
!arrays
 real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size,ndij)
 type(atomorb_type),intent(in) :: Atm
 type(pawrad_type),intent(in) :: Atmrad
 type(pawrad_type),intent(in) :: Pawrad
 type(pawtab_type),target,intent(in) :: Pawtab

!Local variables ---------------------------------------
!scalars
 integer :: ilnc,ilc,lnc_size,l_max
 integer :: lm2_size,ln_size,ln2_size,lmn_size
 integer :: lm_size,klmn,kln,klm
 integer :: lc_max,ilsl,isgg,israd,opt_l,pawprtvol
 real(dp) :: tmp,sgg,intgrl
!character(len=500) :: msg
!arrays
 integer :: opt_l_index(0,0),pack2ij(0)
 integer,allocatable :: kln2ln(:,:),klm2lm(:,:)
 integer, pointer :: indklmn(:,:),indlmn(:,:)
 type(slatang_cshell_t),allocatable :: Slatang3l(:)
 type(slatrad_cshell_t),allocatable :: Slatrad3l(:)

! *************************************************************************

 ! * Consistency check
 ABI_CHECK(cplex_dij==1,"cplex_dij must be 1")

 ABI_CHECK(ndij==1,"ndij must be 1")

 ABI_CHECK(lmn2_size==Pawtab%lmn2_size,"Wrong lmn2_size")

 lmn_size  = Pawtab%lmn_size
 ln_size   = Pawtab%basis_size
 ln2_size  = Pawtab%ij_size
 l_max     = (Pawtab%l_size-1)/2 +1
 lm_size   = l_max**2
 lm2_size  = lm_size*(lm_size+1)/2

 indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
 indklmn => Pawtab%indklmn(1:8,1:lmn2_size)

 ! * Setup of useful tables.
 ABI_MALLOC(kln2ln,(6,ln2_size))
 call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)

 ABI_MALLOC(klm2lm,(6,lm2_size))
 call make_klm2lm(lmn_size,lmn2_size,lm2_size,indlmn,indklmn,klm2lm)

 ! * Integrate angular part.
 lnc_size  = Atm%ln_size
 lc_max    = Atm%l_max

 ABI_DATATYPE_ALLOCATE(Slatang3l,(lm2_size))
 call slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)

 ABI_FREE(klm2lm)

 ! * Integrate radial part.
 ABI_DATATYPE_ALLOCATE(Slatrad3l,(ln2_size))

 call slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad)

 ! === Calculate matrix elements of Sigma_x^core ===
 ! * $<\phi_i|\Sigma_x^\core|\phi_j>$
 dijexc_core = zero
 do klmn=1,lmn2_size
   klm = Pawtab%indklmn(1,klmn)
   kln = Pawtab%indklmn(2,klmn)
   !
   ! * Summing over (lc,nc) and lslat
   tmp = zero
   if (Slatang3l(klm)%nsggsel >0) then
     do ilnc=1,Atm%ln_size
       ilc = 1+Atm%indln(1,ilnc)
       do ilsl=1,Slatang3l(klm)%lslat_max !FIXME check this
         !do ilsl=Slatang3l(klm)%lslat_min,Slatang3l(klm)%lslat_max
         isgg  = Slatang3l(klm)%sggselect(ilsl,ilc)
         israd = Slatrad3l(kln)%rlphic_select(ilsl,ilnc)
         if (isgg>0 .and. israd>0) then
           sgg    = Slatang3l(klm)%sgg(isgg)
           intgrl = Slatrad3l(kln)%rlphic_int(israd)
           tmp = tmp + intgrl*sgg
         end if
       end do
     end do
   end if

   dijexc_core(klmn,1) = -tmp  ! Store results.
 end do

 ! * Print values
 call wrtout(std_out,"   ************** Dij Fock_core ************ ",'COLL')
 opt_l=-1; pawprtvol=-1
 call pawio_print_ij(std_out,dijexc_core(:,1),lmn2_size,cplex_dij,lmn_size,opt_l,opt_l_index,0,pawprtvol,pack2ij,-one,1)

 ! * Free memory.
 ABI_FREE(kln2ln)
 call slatang_cshell_free(Slatang3l)
 ABI_DT_FREE(Slatang3l)
 call slatrad_cshell_free(Slatrad3l)
 ABI_DT_FREE(Slatrad3l)

end subroutine paw_sigxcore
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/paw_mkdijexc_core
!! NAME
!!  paw_mkdijexc_core
!!
!! FUNCTION
!!  Driver routine to calculate the onsite matrix element of the Fock operator between two
!!  all-electron partial waves.
!!
!! INPUTS
!!  ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
!!  cplex_dij=1 if dijexc_core is real, 2 if they are complex
!!  lmn2_size_max=Max Number of (klmn) channels over type of atoms.
!!  Cryst<crystal_t>=Structure describing the crystal structure and its symmmetries.
!!  Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
!!  Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
!!  pawprtvol=Flags governing the verbosity of the output.
!!  filpsp(ntypat)=names of the files containing the all-electron core WF
!!
!! OUTPUT
!!  dijexc_core(cplex_dij*lmn2_size_max,ndij,ntypat)= On-site matrix elements $ \<\phi_i|Sigma_x^\core|\phi_j\>
!!    for each type of atom.
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,pawprtvol,filpsp)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max
 type(crystal_t),intent(in) :: Cryst
!arrays
 real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat) !TODO use ragged arrays pawij?
 character(len=fnlen) :: filpsp(Cryst%ntypat)
 type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)

!Local variables ---------------------------------------
!scalars
 integer :: itypat,ic,ierr,lmn_size,lmn2_size,ln_size,isppol
 real(dp) :: rcut
 character(len=500) :: header,msg
 character(len=fnlen) :: fcore,string
!arrays
 integer,allocatable :: phi_indln(:,:)
 real(dp),ABI_CONTIGUOUS pointer :: phi(:,:)
 real(dp),allocatable :: overlap(:,:)
 type(atomorb_type),allocatable :: Atm(:)
 type(Pawrad_type),allocatable :: Radatm(:)

! *************************************************************************

 ABI_DATATYPE_ALLOCATE(Atm,(Cryst%ntypat))
 ABI_DATATYPE_ALLOCATE(Radatm,(Cryst%ntypat))

 ABI_CHECK(ndij==1     ,"spinor+HF not available")
 ABI_CHECK(cplex_dij==1,"spinor+HF not available")
 ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")

 !allocate(dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) !TODO use ragged arrays pawij?
 dijexc_core=zero

 do itypat=1,Cryst%ntypat

   ! Read core orbitals for this atom type.
   string = filpsp(itypat)
   fcore = "CORE_"//TRIM(basename(string))
   ic = INDEX (TRIM(string), "/" , back=.TRUE.) ! if string is a path, prepend path to fcore.
   if (ic>0 .and. ic<LEN_TRIM(string)) fcore = filpsp(itypat)(1:ic)//TRIM(fcore)

   rcut=Pawtab(itypat)%rpaw
   call init_atomorb(Atm(itypat),Radatm(itypat),rcut,fcore,pawprtvol,ierr)

   if (ierr/=0) then
     msg = " Error reading core orbitals from file: "//TRIM(fcore)
     MSG_ERROR(msg)
   end if
   write(header,'(a,i4,a)')" === Atom type = ",itypat," === "
   call print_atomorb(Atm(itypat),header,unit=std_out,prtvol=pawprtvol)
   !
   ! * Calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for this atom type.
   lmn_size  = Pawtab(itypat)%lmn_size
   lmn2_size = Pawtab(itypat)%lmn2_size

   call paw_sigxcore(cplex_dij,lmn2_size,ndij,&
&    Pawrad(itypat),Pawtab(itypat),Atm(itypat),Radatm(itypat),dijexc_core(1:lmn2_size,:,itypat))

   ln_size =  Pawtab(itypat)%basis_size
   phi     => Pawtab(itypat)%phi

   ABI_MALLOC(phi_indln,(2,ln_size))
   call make_indln(lmn_size,ln_size,Pawtab(itypat)%indlmn(:,:),phi_indln)

   ABI_MALLOC(overlap,(Atm(itypat)%ln_size,ln_size))
   isppol=1 ! hardcoded
   call get_overlap(Atm(itypat),Radatm(itypat),Pawrad(itypat),isppol,ln_size,phi,phi_indln,overlap)

   ABI_FREE(phi_indln)
   ABI_FREE(overlap)
 end do ! ntypat

 ! Free memory
 call pawrad_free(Radatm)
 do itypat=1,Cryst%ntypat
   call destroy_atomorb(Atm(itypat))
 end do

 ABI_DATATYPE_DEALLOCATE(Atm)
 ABI_DATATYPE_DEALLOCATE(Radatm)

end subroutine paw_mkdijexc_core
!!****

!----------------------------------------------------------------------

!!****f* m_paw_slater/slatrad_free_0D
!! NAME
!!  slatrad_free_0D
!!
!! FUNCTION
!!  Free the dynamic memory allocated in a structure of type slatrad_t
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatrad_free_0D(Slatrad)

 implicit none

!Arguments ------------------------------------
!scalars
 type(slatrad_t),intent(inout) :: Slatrad

! *********************************************************************

 !@slatrad_t
 if (allocated(Slatrad%intgl_select)) then
   ABI_FREE(Slatrad%intgl_select)
 end if
 if (allocated(Slatrad%intgl)) then
   ABI_FREE(Slatrad%intgl)
 end if

end subroutine slatrad_free_0D
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/slatrad_free_1D
!! NAME
!!  slatrad_free_1D
!!
!! FUNCTION
!!  Free the dynamic memory allocated in a structure of type slatrad_t
!!
!! PARENTS
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatrad_free_1D(Slatrad)

 implicit none

!Arguments ------------------------------------
!scalars
 type(slatrad_t),intent(inout) :: Slatrad(:)

!Local variables-------------------------------
 integer :: ii
! *********************************************************************

 do ii=1,SIZE(Slatrad)
   call slatrad_free_0D(Slatrad(ii))
 end do

end subroutine slatrad_free_1D
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/slatrad_init
!! NAME
!!  slatrad_init
!!
!! FUNCTION
!!  Initialize the structure storing the radial part of Slater"s integrals.
!!
!! INPUTS
!!  which_intg= Option defining what kind of integrals have to be calculated:
!!   -- 1 for Slater integral of AE partial waves only.
!!          A = \frac{4\pi}{2L+1} \int u_i(1) u_j(1) \frac{r_<^L}{r_^{L+1}} u_k(2)u_l(2) d1d2
!1
!!   -- 2 for Slater integral of (AE-PS) partial waves
!!          B = \frac{4\pi}{2L+1} \int  u_i(1)  u_j(1) \frac{r_<^L}{r_^{L+1}}  u_k(2)  u_l(2) d1d2 -
!!              \frac{4\pi}{2L+1} \int tu_i(1) tu_j(1) \frac{r_<^L}{r_^{L+1}} tu_k(2) tu_l(2) d1d2
!!
!!   -- 3 for Slater integral of (AE-PS-compensation charges)
!!          C = A -
!!          \frac{4\pi}{2L+1} \int [tu_i(1) tu_j(1) + qhat^L_\ij r_1^2 g^L(1) ] \frac{r_<^L}{r_^{L+1}}*
!!                                 [tu_k(2) tu_l(2) + qhat^L_\kl r_2^2 g^L(2) ] d1d2
!!
!!  where u = \phi/r; tu = \tphi/r; and qhat^L_\ij are related to q^\LM_\ij via
!!    q^\LM_\ij = \Gaunt_\ij^\LM qhat^L_\ij => qhat^L\ij = \int (u_i*u_j - tu_i*tu_j) r^L
!!
!!  ln2_size=Number of symmetrical (l,n) channels for this atom type type.
!!  Pawrad<pawrad_type>=paw radial mesh and related data
!!  Pawtab<pawtab_type>=paw tabulated starting data
!!
!! OUTPUT
!!  Slatrad4<slatrad_t>=The object completely initialized.
!!
!! NOTES
!!  Slater integrals S_ij are invariant under exchage of the indeces,
!!  but the results reported by calc_slatradl are not due to numerical roundoff errors (err < 10^-9).
!!  However this does not cause any problem since only the upper triangle of the S_ij matrix
!!  is stored and used in the other routines.
!!
!! PARENTS
!!      m_paw_slater
!!
!! CHILDREN
!!      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout
!!
!! SOURCE

subroutine slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad,Pawtab)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: ln2_size,which_intg
!arrays
 type(pawrad_type),target,intent(in) :: Pawrad
 type(pawtab_type),target,intent(in) :: Pawtab
 type(slatrad_t),intent(out) :: Slatrad4(ln2_size*(ln2_size+1)/2)

!Local variables ---------------------------------------
!scalars
 integer :: mesh_size,il,iln,isl,in,jl,jln,jn,sln1,sln2,l_slat
 integer :: kn,kl,ln,ll,kln,lln,lslat_max,lslat_min,nintgl
 integer :: lmn_size,lmn2_size,ln_size,slt_idx
 real(dp) :: ae_intg,ps_intg,pshat_intg,intg,tqij_L,tqkl_L !intg1
 character(len=500) :: msg
!arrays
 integer,allocatable :: kln2ln(:,:)
 integer, pointer :: indklmn(:,:),indlmn(:,:)
 real(dp),allocatable :: uiuj(:),ukul(:),tuituj(:),tuktul(:),tuituj_tqgl(:),tuktul_tqgl(:)
 real(dp),allocatable :: tmp_integrals(:),ff(:)
 real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:),phi_k(:),phi_l(:)
 real(dp),ABI_CONTIGUOUS pointer :: tphi_i(:),tphi_j(:),tphi_k(:),tphi_l(:)
 real(dp),ABI_CONTIGUOUS pointer :: shapefunc(:),rad(:)

! *************************************************************************

 DBG_ENTER("COLL")

 ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")

 if ( ALL(which_intg /= (/1,2,3/)) ) then
   write(msg,'(a,i0)')"Wrong value for which_intg: ",which_intg
   MSG_ERROR(msg)
 end if

 !@slatrad_t
 lmn_size   = Pawtab%lmn_size
 lmn2_size  = Pawtab%lmn2_size
 ln_size    = Pawtab%basis_size
 mesh_size  = Pawtab%mesh_size
 !
 ! Useful table for looping.
 indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
 indklmn => Pawtab%indklmn(1:8,1:lmn2_size)

 ABI_MALLOC(kln2ln,(6,ln2_size))
 call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)

 ABI_MALLOC(uiuj,(mesh_size))
 ABI_MALLOC(ukul,(mesh_size))
 ABI_MALLOC(ff,(mesh_size))
 ABI_MALLOC(tuituj,(mesh_size))
 ABI_MALLOC(tuktul,(mesh_size))
 ABI_MALLOC(tuituj_tqgl,(mesh_size))
 ABI_MALLOC(tuktul_tqgl,(mesh_size))
 rad => Pawrad%rad
 !
 ! * Loop over (k,l) channels in packed form.
 do sln2=1,ln2_size
   kl  = kln2ln(1,sln2)
   ll  = kln2ln(2,sln2)
   kn  = kln2ln(3,sln2)
   ln  = kln2ln(4,sln2)
   kln = kln2ln(5,sln2)
   lln = kln2ln(6,sln2)
   !write(std_out,*)"sln2, kln, lln",sln2,kln,lln

   phi_k  => Pawtab%phi (:,kln)
   tphi_k => Pawtab%tphi(:,kln)

   phi_l  => Pawtab%phi (:,lln)
   tphi_l => Pawtab%tphi(:,lln)
   !
   ! * Loop over (i,j) channels in packed form AND only for the upper triangle of sln2, sln1
   do sln1=1,sln2
     il  = kln2ln(1,sln1)
     jl  = kln2ln(2,sln1)
     in  = kln2ln(3,sln1)
     jn  = kln2ln(4,sln1)
     iln = kln2ln(5,sln1)
     jln = kln2ln(6,sln1)
     !write(std_out,*)"sln1, iln, jln",sln1,iln,jln

     slt_idx = sln1 + sln2*(sln2-1)/2 ! index for packed storage.

     phi_i  => Pawtab%phi (:,iln)
     tphi_i => Pawtab%tphi(:,iln)

     phi_j  => Pawtab%phi (:,jln)
     tphi_j => Pawtab%tphi(:,jln)

     lslat_min = MAX(ABS(il-jl),ABS(kl-ll)) + 1  ! We use indeces not l-values.
     lslat_max = MIN((il+jl),(kl+ll)) - 1

     !lslat_min = MIN(ABS(il-jl),ABS(kl-ll)) + 1
     !lslat_max = MAX((il+jl),(kl+ll)) - 1

     Slatrad4(slt_idx)%lslat_min = lslat_min
     Slatrad4(slt_idx)%lslat_max = lslat_max

     Slatrad4(slt_idx)%iln = iln
     Slatrad4(slt_idx)%jln = jln
     Slatrad4(slt_idx)%kln = kln
     Slatrad4(slt_idx)%lln = lln

     ABI_MALLOC(Slatrad4(slt_idx)%intgl_select,(lslat_min:lslat_max))
     Slatrad4(slt_idx)%intgl_select(:) = 0
     Slatrad4(slt_idx)%nintgl          = 0

     if (lslat_min > lslat_max) then
       ! e.g. (1 2) (1 1). Due to angular selection rules, this integral do not contribue
       !write(std_out,*)"lslat_min, lslat_max",lslat_min,lslat_max
       !write(std_out,*)"il,jl,kl,ll",il,jl,kl,ll
       !MSG_ERROR("")
       ABI_MALLOC(Slatrad4(slt_idx)%intgl,(0))
       CYCLE
     end if

     uiuj   =  phi_i *  phi_j  ! The AE part.
     ukul   =  phi_k *  phi_l
     tuituj = tphi_i * tphi_j  ! The pseudized part.
     tuktul = tphi_k * tphi_l
     !
     ! Calculate L-depedent integrals where L come from the expansion the Coulomb interaction.
     ABI_MALLOC(tmp_integrals,(MAX(lslat_min,lslat_max)))
     tmp_integrals=zero
     nintgl=0

     do isl=lslat_min,lslat_max
     !do isl=lslat_min,lslat_max,2  ! TODO Here I can reduce the number of iterations using a step of 2.
       l_slat = isl-1
       call calc_slatradl(l_slat,mesh_size,uiuj,ukul,Pawrad,ae_intg)
       intg = ae_intg

#if 0
       call calc_slatradl(l_slat,mesh_size,ukul,uiuj,Pawrad,intg1)
       if (ABS(intg1-ae_intg)>tol12) then
         write(msg,'(a,es16.8)')"s_ij and s_ij differ by ",intg1-ae_intg
         MSG_WARNING(msg)
       end if
#endif
       if (which_intg == 2) then ! Subtract the pseudo part.
         call calc_slatradl(l_slat,mesh_size,tuituj,tuktul,Pawrad,ps_intg)
         intg = intg - ps_intg

       else if (which_intg == 3) then ! Subtract (pseudo + compensation charges)
         !
         ! Evaluate tqij_L and tqkl_L (without M-dependent part).
         ff(1)=zero
         ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln1)-pawtab%tphitphj(2:mesh_size,sln1))*rad(2:mesh_size)**l_slat
         if (l_slat==0.and.kl==1.and.ll==1) then
           call pawrad_deducer0(ff,mesh_size,pawrad)
         end if
         call simp_gen(tqij_L,ff,pawrad)

         ff(1)=zero
         ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln2)-pawtab%tphitphj(2:mesh_size,sln2))*rad(2:mesh_size)**l_slat
         if (l_slat==0.and.il==1.and.jl==1) then
           call pawrad_deducer0(ff,mesh_size,pawrad)
         end if
         call simp_gen(tqkl_L,ff,pawrad)

         shapefunc   => Pawtab%shapefunc(:,isl)  ! Recheck this part, in particular the convention
         tuituj_tqgl = tuituj + tqij_L * shapefunc * rad**2
         tuktul_tqgl = tuktul + tqkl_L * shapefunc * rad**2

         call calc_slatradl(l_slat,mesh_size,tuituj_tqgl,tuktul_tqgl,Pawrad,pshat_intg)
         intg = intg - pshat_intg
       end if
       !
       ! * Store results
       if (ABS(intg)>=tol12) then
         nintgl = nintgl +1
         Slatrad4(slt_idx)%intgl_select(isl) = nintgl
         tmp_integrals(nintgl) = intg
       end if
     end do !isl
     !
     ! Finalize the object.
     Slatrad4(slt_idx)%nintgl = nintgl
     ABI_MALLOC(Slatrad4(slt_idx)%intgl,(nintgl))
     if (nintgl>0) Slatrad4(slt_idx)%intgl(:) = tmp_integrals(1:nintgl)
     ABI_FREE(tmp_integrals)
   end do !sln1
 end do !sln2
 !
 ! Free memory
 ABI_FREE(kln2ln)
 ABI_FREE(uiuj)
 ABI_FREE(ukul)
 ABI_FREE(tuituj)
 ABI_FREE(tuktul)
 ABI_FREE(ff)
 ABI_FREE(tuituj_tqgl)
 ABI_FREE(tuktul_tqgl)

 DBG_EXIT("COLL")

end subroutine slatrad_init
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/paw_dijhf
!! NAME
!!  paw_dihf
!!
!! FUNCTION
!!  This routine calculates the onsite D_{ij} strengths of the exchange parth of the self energy.
!!
!! INPUTS
!!  ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
!!  cplex_dij=1 if sigx_dij is real, 2 if they are complex
!!  qphase=2 if dij contains a exp(-i.q.r) phase (as in the q<>0 RF case), 1 if not
!!  lmn2_size_max=Max Number of (klmn) channels over type of atoms.
!!  my_natom=number of atoms treated by current process
!!  ntypat=number of atom types
!!  Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
!!  Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
!!  Pawang<type(pawang_type)>=paw angular mesh and related data
!!  pawprtvol=Flags governing the verbosity of the output.
!!
!! OUTPUT
!!  sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom)=
!!    For each atom, the Pseudopotential strengths of the on-site operator Sigma_x
!!
!! NOTES
!!  The on-site contribution to the matrix elements of the exchange part of the self-energy is given by:
!!  <\tpsi_a| [\sum_{ij} |tprj_i\> D_{ij} \<tprj_j|] |\tpsi_b\>.
!!
!!  When compensation charges are used one obtains:
!!
!!  D_{ij} = - sum_{kl} \rho_lk [ \Phi_{ikjl} - \Phihat_{ijkl} =
!!
!!         = - sum_{kl} \rho_lk \sum_{LM} \Gaunt_{ik}^{LM} \Gaunt_{jl}^{LM} [S_{ikjl}^L - tS{ikjl}^L}]
!!
!!  where S and tS are Slater-like integrals given by
!!
!!  1)  S_{ijkl}^L = dfrac{4\pi}{2L+1} \iint u_i(1)u_j(1) u_k(2) u_l(2) \dfrac{r_<^L/}{r_>^{L+1}} d1d2.
!!  1) tS_{ijkl}^L = dfrac{4\pi}{2L+1} \iint [u_i(1)u_j(1)+ tq_{ij}^L g^L(1)]
!!                                           [u_k(2)u_l(2)+ tq_{kl}^L g^L(2)] \dfrac{r_<^L/}{r_>^{L+1}} d1d2.
!!
!!  tq_{ij}^L is defined in terms of q_{ij}^L via: q_{ij]^{LM} = tq_{ij}^L \Gaunt_{ij}^{LM}
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

subroutine paw_dijhf(ndij,cplex_dij,qphase,lmn2_size_max,my_natom,ntypat,Pawtab,Pawrad,Pawang,Pawrhoij,&
&                    sigx_dij,pawprtvol)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max,my_natom,ntypat,qphase
 type(pawang_type),intent(in) :: Pawang
!arrays
 real(dp),target,intent(out) :: sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom) !TODO use ragged arrays pawij?
 type(pawtab_type),intent(in) :: Pawtab(ntypat)
 type(pawrad_type),intent(in) :: Pawrad(ntypat)
 type(pawrhoij_type),intent(in) :: Pawrhoij(my_natom)

!Local variables ---------------------------------------
!scalars
 integer,parameter :: cplex=1    ! FIXME preliminary implementation
 integer :: cplex_rhoij,iatom,iq,iq0_dij,iq0_rhoij,itypat,lmn_size,lmn2_size,ispden,nspden,ln2_size
 integer :: lm2_size !,isppol ln_size,
 integer :: irhoij,jrhoij
 integer :: rho_lmn !,rho_klm,rho_kln,rho_lmin,rho_lmax,rho_iln,rho_jln
 integer :: klmn
 integer :: i_lmn,j_lmn,k_lmn,l_lmn
 integer :: which_intg,l_max,opt_l
 real(dp) :: ro,slt_ikjl,slt_iljk
 !character(len=500) :: msg
!arrays
 integer :: opt_l_index(0,0),pack2ij(0)
 real(dp), ABI_CONTIGUOUS pointer :: sigx_atm(:,:)
 type(slatrad_t),allocatable :: Slatrad4(:)

! *************************************************************************

 DBG_ENTER("COLL")

 ABI_CHECK(ndij/=4,"ndij=4 not coded")
 ABI_CHECK(cplex_dij==1,"cplex_dij/=2 not coded")
 ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")

 if (my_natom>0) then
   if (pawrhoij(1)%qphase<cplex) then
     MSG_BUG('Must have pawrhoij()%qphase >= cplex !')
   end if
 end if

 sigx_dij=zero

 do iatom=1,my_natom
   itypat   =Pawrhoij(iatom)%itypat
   lmn_size =Pawtab(itypat)%lmn_size
   lmn2_size=Pawtab(itypat)%lmn2_size
   l_max    =(Pawtab(itypat)%l_size+1)/2
   lm2_size = (l_max**2)*(l_max**2+1)/2
   cplex_rhoij=Pawrhoij(iatom)%cplex_rhoij
   !write(std_out,*)"in atom ",iatom,"lm2_size=",lm2_size

   ! Calculate Slater integral for this atom type.
   ! TODO obviously these tables should be stored in Pawtab!
   ln2_size = Pawtab(itypat)%ij_size
   ABI_DATATYPE_ALLOCATE(Slatrad4,(ln2_size*(ln2_size+1)/2))
   which_intg=3
   call slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad(itypat),Pawtab(itypat))

   sigx_atm => sigx_dij(:,:,iatom)

!  Loop over phase exp(iqr) phase real/imaginary part, if any
   do iq=1,qphase
     !First loop: we store the real part in dij(1 -> lmn2_size)
     !2nd loop: we store the imaginary part in dij(lmn2_size+1 -> 2*lmn2_size)
     iq0_dij=merge(0,cplex_dij*lmn2_size,iq==1)
     iq0_rhoij=merge(0,cplex_rhoij*lmn2_size,iq==1)

     ! * Loop over spin components.
     nspden=ndij
     do ispden=1,ndij
       !
       ! ============================================================
       ! ==== Summing over the non-zero lk channels of \rho_{lk} ====
       ! ============================================================
       jrhoij=1+iq0_rhoij
       do irhoij=1,pawrhoij(iatom)%nrhoijsel
         rho_lmn=pawrhoij(iatom)%rhoijselect(irhoij)

         ! check wheter rho_lmin is consistent with the Indexing used in slatrad
         !rho_klm =pawtab(itypat)%indklmn(1,rho_lmn)
         !rho_kln =pawtab(itypat)%indklmn(2,rho_lmn)
         !rho_lmin=pawtab(itypat)%indklmn(3,rho_lmn)
         !rho_lmax=pawtab(itypat)%indklmn(4,rho_lmn)

         ! Retrieve rhoij for this ispden.
         if (nspden/=2) then
           ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
         else
           MSG_ERROR("Recheck this part")
           if (ispden==1) then
             ro=pawrhoij(iatom)%rhoijp(jrhoij,1) + pawrhoij(iatom)%rhoijp(jrhoij,2)
           else if (ispden==2) then
             ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
           end if
         end if
         !
         ! Avoid double-counting the diagonal of rho.
         ro=ro*pawtab(itypat)%dltij(rho_lmn)*half

         call klmn2ijlmn(rho_lmn,lmn_size,k_lmn,l_lmn)

         ! Loop over the upper triangle of the D_{ij) matrix and accumulate:
         ! sum_\lk rho_\kl [ \Phi_{ikjl} + \Phi_{iljk} - \Phihat_{ikjl} - \Phihat_{iljk} ]
         do klmn=1,lmn2_size
           ! Calculate the indeces in the Slatrad4 structure.
           call klmn2ijlmn(klmn,lmn_size,i_lmn,j_lmn)

           ! My formula
           slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,j_lmn,l_lmn)
           slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,j_lmn,k_lmn)

           !slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,l_lmn,j_lmn)
           !slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,k_lmn,j_lmn)
           !slt_iljk = slt_ikjl

           sigx_atm(klmn+iq0_dij,ispden) = sigx_atm(klmn,ispden) + ro * (slt_ikjl + slt_iljk)
         end do ! klmn

         jrhoij=jrhoij+cplex_rhoij
       end do ! irhoij
     end do ! iq
   end do ! ispden

   if (ABS(pawprtvol)>=1) then ! * Print values
     call wrtout(std_out,"   ************** Dij Fock ************ ",'COLL')
     opt_l=-1
     call pawio_print_ij(std_out,sigx_atm(:,1),lmn2_size,cplex_dij,lmn_size,opt_l,opt_l_index,0,pawprtvol,pack2ij,-one,1)
   end if

   call slatrad_free(Slatrad4)
   ABI_DATATYPE_DEALLOCATE(Slatrad4)
 end do ! iatom

 ! Factor half cancels in the derivation wrt rho_ij.
 sigx_dij = - sigx_dij

 DBG_EXIT("COLL")

end subroutine paw_dijhf
!!****

!----------------------------------------------------------------------

!!****f* m_paw_slater/summ_2gaunt
!! NAME
!!   summ_2gaunt
!!
!! FUNCTION
!!  Helper function returning \sum_M G_{ij}^{LM} G_{kl}^{LM}
!!
!! INPUTS
!!  ij_lm=index of (i_lm,j_lm) element in packed form.
!!  kl_lm=index of (k_lm,l_lm) element in packed form.
!!  ll_idx=Index for L (thus L+1).
!!  Pawang<type(pawang_type)>=paw angular mesh and related data
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

function summ_2gaunt(Pawang,ij_lm,kl_lm,ll_idx)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: ij_lm,kl_lm,ll_idx
 real(dp) :: summ_2gaunt
 type(pawang_type),intent(in) :: Pawang
!arrays

!Local variables-------------------------------
!scalars
 integer :: ignt1,ignt2,idx_LM,max_klm,mm,ii,ll
 character(len=500) :: msg

!************************************************************************

 ! FIXME: size of gntselect depends on pawxcdev!
 ! Consistency check on input arguments.
 max_klm = pawang%l_max**2*(pawang%l_max**2+1)/2
 if (ij_lm>max_klm.or.kl_lm>max_klm.or.ij_lm<1.or.kl_lm<1.or.&
&    ll_idx>pawang%l_size_max.or.ll_idx<1) then
   write(msg,'(a,3i0)')"Wrong indeces, check pawxcdev ",ij_lm,kl_lm,ll_idx
   MSG_ERROR(msg)
 end if

 ll = ll_idx-1
 summ_2gaunt=zero; ii=0
 do mm=-ll,ll
   idx_LM = 1 + ll**2 + ll + mm
   ignt1 = Pawang%gntselect(idx_LM,ij_lm)
   ignt2 = Pawang%gntselect(idx_LM,kl_lm)
   if (ignt1>0 .and. ignt2>0) then
     summ_2gaunt = summ_2gaunt + Pawang%realgnt(ignt1)*Pawang%realgnt(ignt2)
     ii=ii+1
     write(std_out,'(a,4(i2,1x),f8.5,i2)')"ll, mm, ij_lm, kl_lm: ",ll,mm,ij_lm,kl_lm,summ_2gaunt,ii
     if (ii/=1) MSG_WARNING("ii>1")
   end if
 end do

end function summ_2gaunt
!!***

!----------------------------------------------------------------------

!!****f* m_paw_slater/slat_intg
!! NAME
!!   slat_intg
!!
!! FUNCTION
!!  Helper function returning the slater integral
!!    \int_\Omega \phi_i(1)\phi_j(1) \dfrac{1}{|1-2|} \phi_k(2)\phi_l(2) d1d2
!!
!! INPUTS
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

function slat_intg(Slatrad4,Pawtab,Pawang,i_lmn,j_lmn,k_lmn,l_lmn)

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: i_lmn,j_lmn,k_lmn,l_lmn
 real(dp) :: slat_intg
 type(pawtab_type),intent(in) :: Pawtab
 type(pawang_type),intent(in) :: Pawang
!arrays
 type(slatrad_t),intent(in) :: Slatrad4(:)

!Local variables-------------------------------
!scalars
 integer :: ij_lmn,kl_lmn,kl_ln,ij_lm,kl_lm,ilsum,ij_ln
 integer :: isel,slt_idx
 integer :: iln,jln,kln,lln,ii
 real(dp) :: sltL_ijkl,angintL_ijkl
 !character(len=500) :: msg

!************************************************************************

 ! The lmn packed indeces for (ij) and (kl).
 if (j_lmn>=i_lmn) then
   ij_lmn = i_lmn + j_lmn*(j_lmn-1)/2
 else
   ij_lmn = j_lmn + i_lmn*(i_lmn-1)/2
 end if

 if (l_lmn>=k_lmn) then
   kl_lmn = k_lmn + l_lmn*(l_lmn-1)/2
 else
   kl_lmn = l_lmn + k_lmn*(k_lmn-1)/2
 end if
 !
 ! The lm indeces for (ij) and (kl) in packed storage.
 ij_lm = pawtab%indklmn(1,ij_lmn)
 ij_ln = pawtab%indklmn(2,ij_lmn)

 kl_lm = pawtab%indklmn(1,kl_lmn)
 kl_ln = pawtab%indklmn(2,kl_lmn)
 !
 ! The index of (ijkl) in the Slatrad4 database.
 if (kl_ln>=ij_ln) then
   slt_idx = ij_ln +kl_ln*(kl_ln-1)/2
 else
   slt_idx = kl_ln +ij_ln*(ij_ln-1)/2
 end if

!BEGIN DEBUG
 iln = Slatrad4(slt_idx)%iln
 jln = Slatrad4(slt_idx)%jln
 kln = Slatrad4(slt_idx)%kln
 lln = Slatrad4(slt_idx)%lln

 ii = kln + lln*(lln-1)/2
 if (slt_idx /=  (iln + jln*(jln-1)/2 + ii*(ii-1)/2 )) then
   write(std_out,*)"slt_idx, iln, jln, kln, lln",slt_idx, iln, jln, kln, lln
   MSG_BUG("Check indeces")
 end if
!END DEBUG
 !
 ! Calculate the integral by summing over ilsum.
 slat_intg=zero
 if (Slatrad4(slt_idx)%nintgl>0) then
   do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max
   !% do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
     isel = Slatrad4(slt_idx)%intgl_select(ilsum)
     if (isel/=0) then
       sltL_ijkl = Slatrad4(slt_idx)%intgl(isel)
       angintL_ijkl = summ_2gaunt(Pawang,ij_lm,kl_lm,ilsum)
       slat_intg = slat_intg + sltL_ijkl * angintL_ijkl
     end if
   end do
 end if

end function slat_intg
!!***

!----------------------------------------------------------------------

END MODULE m_paw_slater
!!***
